shifted QR algorithm implementation for computing roots given coefficients by pasquale90 · Pull Request #3 · Audio-Project-Workgroup/digital_filter_visualiser

pasquale90 · 2026-06-10T12:03:49Z

@acknowledgment : Many thanks to @LizAryslanova and my sincere gratitude for her work in analyzing the algorithm and for her guidance in helping me understand it

Method used to convert polynomial coefficients into roots using a QR method.
Consists of the following steps:

Preprocessing:
Filter out leading zeros if any
Filter out trailing zeros if any
Build companion matrix of Hessenberg form from given coefficients
Calculate rayleigh quotient shift
apply shifted QR iterations (starting from bottom right) with deflation when subdiagonal entries become ~0
- Step 1 calculate Q0 via Gram Schmidt orthogonalization method
- Step 2 calculate R0 : R0 = Q0 A0
- Step 3 calculate A1 : A1 = R0 Q0
- Step 4 add shift back in A1 and check sub-diagonal entries for deflation
Extract eigenvalues:
- merge roots that have smaller scaled difference than tolerance
- scan subdiagonal:
  - if values on the subdiagonal smaller than Epsilon, then real root
  - else 2 solutions. Solve det(A - λI) = 0 for 2x2 square matrix:
    - if discriminant >= 0, then two real solutions
    - else then root has a complex conjugate (conjugate is not appended on returned vector, as this is done automatically when using FilterState::add to add that root in the Value tree)
Returns complex roots paired with their corresponding order.

Leftover tasks:

Fix bug when editing intermediate zero-valued coefficients in between leading or trailing zero-valued coefficients.
Implement tests
- For assessing level of accuracy (distance and order)
- For assessing performance (execution time)
Evaluate other shifting techniques (double shift, wilkinson shift)

justonem0reuser

Thank you,
All the existing tests are passed.
Unfortunately, I didn't go into details of the method itself, so I can only propose some ways for optimization.

Have you considered using 1D vectors for A, Q, R instead of "vectors of vectors"?
I believe this would speed up the process significantly.
The drawback is the necessity of manual indexing, for example: Q[i * degree + j] instead of Q[i][j].
I wouldn't perform memory allocation in cycle.
The "old-school" optimization (though probably it doesn't make sense nowadays): not to repeat some calculations (sqrt) and prefer multiplication to division.
Matrix multiplication can also be accelerated, but it requires more sophisticated optimization.
In TODO note it is proposed to change QR return type to bool for the case when there is a pole outside the unit. Is it really necessary? CoefficientsToRoots is an abstract math class that doesn't solve only the particular filter design task, so it is probably not its responsibility to perform this check.

justonem0reuser · 2026-06-10T22:59:01Z

+        for (size_t col = 0 ; col < shift_idx; col++)
+        {
+            // set v with column A[col]
+            std::vector<double> v(shift_idx);


Suggested change

std::vector<double> v(shift_idx);

OK, the allocation of this vector is now moved outside the loops e75f2ae

justonem0reuser · 2026-06-10T23:01:06Z

+    std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
+    std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
+
+    size_t shift_idx {degree}, iter {0};


Suggested change

size_t shift_idx {degree}, iter {0};

size_t shift_idx {degree}, iter {0};

std::vector<double> v(shift_idx);

justonem0reuser · 2026-06-10T23:04:05Z

+            norm = std::sqrt(norm);
+
+            for (size_t k = 0; k < shift_idx; ++k)
+            {
+                Q[k][col] = v[k] / norm;
+            }


Suggested change

norm = std::sqrt(norm);

for (size_t k = 0; k < shift_idx; ++k)

{

Q[k][col] = v[k] / norm;

}

double normReciprocal = 1.0 / std::sqrt(norm);

for (size_t k = 0; k < shift_idx; ++k)

{

Q[k][col] = v[k] * normReciprocal;

}

OK, fixed that here f4bb673

Btw, next time feel free to push the changes yourself , you deserve the credits :)

justonem0reuser · 2026-06-10T23:05:05Z

+            if ( discriminant >= 0.0)
+            {
+                // tow real roots
+                addRoot(c128(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0));
+                addRoot(c128(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0));
+            }
+            else
+            {
+                // Complex conjugate pair
+                const double re = (a+d) / 2.0;
+                const double im = std::sqrt(-discriminant) / 2.0;
+                addRoot(c128(re,im));
+                // addRoot(c128(re,-im)); // Note: this is added automatically later using FilterState::add method. Commenting this, removes the bug of overlapping roots.
+            }


Suggested change

if ( discriminant >= 0.0)

{

// tow real roots

addRoot(c128(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0));

addRoot(c128(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0));

}

else

{

// Complex conjugate pair

const double re = (a+d) / 2.0;

const double im = std::sqrt(-discriminant) / 2.0;

addRoot(c128(re,im));

// addRoot(c128(re,-im)); // Note: this is added automatically later using FilterState::add method. Commenting this, removes the bug of overlapping roots.

}

const double halfSum = 0.5 * (a + d);

const double halfSqrt = 0.5 * std::sqrt(std::abs(discriminant));

if ( discriminant >= 0.0)

{

// tow real roots

addRoot(c128(halfSum + halfSqrt, 0.0));

addRoot(c128(halfSum - halfSqrt, 0.0));

}

else

{

// Complex conjugate pair

const double re = halfSum;

const double im = halfSqrt;

addRoot(c128(re,im));

// addRoot(c128(re,-im)); // Note: this is added automatically later using FilterState::add method. Commenting this, removes the bug of overlapping roots.

}

Fixed as mentioned here #3 (comment)

ry-gatoni · 2026-06-12T21:07:32Z

This is great work Paschalis.

I agree with all of Anton's points and suggestions. Indeed compilers will perform common subexpression elimination in optimized builds, but performing it manually can speed up debug builds and simplify refactoring. The memory layout/ access patterns will be completely reworked during the optimization pass, so there's no point worrying about that for this PR.

The case of a pole outside the unit circle should be handled in the coefficients editor. When you compute the new poles just scan the vector for any unstable poles. If you find any, don't update the filter state. There should be some indication to the user that coefficient they just entered made the pole unstable, though it's hard to give actionable advice for how to supply coefficients which result in stable filters. Currently, trying to add a pole outside the unit circle will not do anything, which means the editor will display different coefficients from what a user enters.

…d on the GramSchmidt method. Computes both real and complex roots Cons: Inefficient as number of coefficients scale up. Optimizations and alternative methods that provide a more efficient way to handle input should be considered for providing a practical application of this feature

…ble and maintainable. + added a todo as a comment for improving the data structure that stores the extracted roots in the upcoming updates

1. fix bug in updateFilterStateOnCoefEdit when adding a pole that already exists (increasing its order), and then removing it (removing the previously added new pole) 2. fix index of returned root when degree equals 1.0 3. fix sign of returned root when degree equals 1.0 4. fix coef_index when constructing the companion matrix. Now it always initializes the matrix with the last degree-lengthed coefficients (correctly ignoring the leading zeros and 1.0)

change termination condition, now all subdiagonal values should be below epsilon (not their sum) fix typo : step 3 - A = R Q change tolerance threshold to 5e-2 add order in prints for debugging

trailing zeros in (feedback) coefficients (when adding default poles) are now parsed out of the companion matrix and roots at 0.,0. (of order equal to #numOfTrailingZeros) are appended directly into the root list before building the companion matrix and applying the QR prev_sz is now calculated directly before appending new poles on the filterstate, taking into account that some updated roots may not result in removing previous poles, but changing their order instead leading zeros if (feedforward) coefficients (when total order of poles exceeds total order of zeros) are now tracked using std::numeric_limits instead of comparing doubles with the == operator.

Fix bug when degree equals to zero. Now the roots (including potential roots at (0,0)), are returned reverse last fix on parsing leading zeros by using std::numeric_limits instead of direct comparison. Now we compare directly using the == operator.

… from GramSchmidt to QR

…). An empirical tolerance variable is used allowing errors of 1 order per 4 orders (scaling per root)

…o-optimization) : Optimization change: replace 2D vectors for A,Q,R with 1D vectors, applying manual indexing ([row][col] now becomes [row * degree + col]). Micro-optimization change : Avoid repeated memory allocation by moving the v vector for storing current columns of A outside loops, and setting its size to {max shift_idx}, that equals to degree of the polynomial

pasquale90 · 2026-06-14T10:39:33Z

@justonem0reuser @ry-gatoni

Thank you for the recommendations . This a very constructive feedback and I totally agree with your points. I was about adding tests before applying any change, but its fine working on them after demystifying the following

regarding tests

Before applying any change, I considered adding tests to have a reference and a basis for comparison. I ve been able to push some tests here (50f60fe). These are comparing the total order of the computed roots to the total order of the given ones, with an empirical tolerance for higher order roots ( 1 order of tolerance for each root, i.e. a 4th order root would add up to the total tolerance by 1, a 8th order root by 2 and so on).

Here are some concerns about the testing framework of this implementation:

TLDR; Single root testing simple and straightforward / multiple root testing black hole and pain.

These tests are now making use of 2 different implementations : Input --> RootsToCoefficients() --> CoefficientsToRoots(). Therefore , as @justonem0reuser suggested, they should be considered integration tests, not unit tests, as they don t test these two functions independently. And I agree. Finding the ground truth, for single roots is straightforward (i.e. single root at {2, -3} has two complex-conjugated roots: (2+3i) and (2-3i) and inversely), but what about when testing multiple roots? Shouldn't we need paper and pen for that? Cause this can become quite cumbersome for extensive tests like these ones. Otherwise do you think these integration tests are sufficient for testing the overall framework (roots -> coeffs and coeffs -> roots)?
Testing the alignment between input roots and computed roots, would require matching roots from one group to roots of the other group and then compute the error for each subgroup to calculate the total error. Unless there is another way (and I really hope it does!). But based on our current approach:

In case of matching computed with actual roots, a quick research led me to the Generalized Assigning Problem or to Bin_packing_problem with some modifications, but requires further searching - and it seems quite challenging.
That s not the case for testing filters with a unique root of arbitrary order. We can have a quite robust and accurate error calculation, or even for a limited amount of roots that are far in the space. Otherwise, a naive solution would be computing the error for each computed root to the nearest given root. Simple but not accurate, especially for given roots that are close to each other, or that has higher order (as computed roots tend to spread more around a root as its order increases).

Metrics:
a) Mean Absolute Error (MAE) could probably be used, but I am not sure if polar coordinates have to be converted to cartesian coordinates before computing the euclidean distance, and still not sure what semantics would that capture.
b) computing (from the angle) the frequency error and phase error, and the exponential error (from the radial difference).
In both cases, we still need to address the case of order misalignment (i.e. input root of order x and computed root of order <x), or when given pole (not a zero) near the unit circle, and computed pole outside of it.

… loop over subdiabonal elements and b) by avoiding double call of sqrt + replacing division with 2.0 by multipying with 0.5 instead on the extractRoots method

pasquale90 · 2026-06-14T11:22:38Z

Have you considered using 1D vectors for A, Q, R instead of "vectors of vectors"?
I believe this would speed up the process significantly.
The drawback is the necessity of manual indexing, for example: Q[i * degree + j] instead of Q[i][j].

Yes, this is a critical optimization step. Replaced 2D vectors for A,Q,R with 1D vectors here e75f2ae

Both prime integration tests and manual tests seem to produce the same results.

…ating causality constraint (at least one pole outside the unit circle - magnitude >= 1.0)

pasquale90 · 2026-06-14T13:20:02Z

The case of a pole outside the unit circle should be handled in the coefficients editor. When you compute the new poles just scan the vector for any unstable poles. If you find any, don't update the filter state. There should be some indication to the user that coefficient they just entered made the pole unstable, though it's hard to give actionable advice for how to supply coefficients which result in stable filters. Currently, trying to add a pole outside the unit circle will not do anything, which means the editor will display different coefficients from what a user enters.

Thanks you @ry-gatoni

In this commit 445fdd6 a check before updating the filterstate implemented, avoiding filterstate updates in case of instability. Informing the user about it, is not yet implemented and left as a TODO item. Lmk if that s properly handled and fine to leave as it is for now

ry-gatoni · 2026-06-15T16:30:30Z

A few notes:

there should be no "order tolerance"; the order is either correct or the test failed.
- note that discrepancies can be handled by the filter state; when you add two roots near enough to each other, they will be combined.
- if the order is still wrong after close-enough roots have been combined, then the method needs to be modified
in light of the above, we can make this a full integration test and check against the filter state roots instead of those returned by the QR method.
there is a bug in TestHelper::makeFilterState; it does not properly clear the filter state. I have a fix for this which I will push to master shortly
test correctness should involve two things:
- do we get the expected number of roots?
  - if not we need to tune the method and/or the filter state merge tolerance
- are the roots in the right place?
  - should use absolute distance, ie std::abs(computed - expected), or maybe relative std::abs(computed - expected) / std::abs(expected). don't worry about root magnitude and phase and how they affect the frequency response; if this method is accurate enough, they will be too.

So I don't think this should be a unit test per se; the metric for "correctness" should be "do the roots end up in the expected place with the expected orders when you add them to the filter state", since that's what we care about at the end of the day anyway.

As for the case of unstable poles, we can discuss how to notify the user on Wednesday. It's fine not having that done for this PR.

Once the tests have been modified to check what I described above and they are passing, I think this PR will be good to merge.

I'm happy to contribute some of the modifications I've made while reviewing these changes if that would be helpful.

justonem0reuser reviewed Jun 10, 2026

View reviewed changes

pasquale90 added 13 commits June 13, 2026 13:15

improve extract roots method, symplifying it and making it more reada…

1241d52

…ble and maintainable. + added a todo as a comment for improving the data structure that stores the extracted roots in the upcoming updates

fix typo in CoefficientsToRoots::extractRoots method

6f78561

implement updating of the filter state on coefficient edits

e1fcfe8

minor improvements and fixes:

f0dc2f2

change termination condition, now all subdiagonal values should be below epsilon (not their sum) fix typo : step 3 - A = R Q change tolerance threshold to 5e-2 add order in prints for debugging

supress gcc warnings

2e9d98b

implement rayleigh quotient shift to speed up QR convergence

6cd4ac6

2 fixes:

307dd97

Fix bug when degree equals to zero. Now the roots (including potential roots at (0,0)), are returned reverse last fix on parsing leading zeros by using std::numeric_limits instead of direct comparison. Now we compare directly using the == operator.

add documentation on CoefficientsToRoots class and change method name…

5332c6d

… from GramSchmidt to QR

clean print statements

b2d6f3a

add tests for Coefficients2Roots::QR function (order computation only…

50f60fe

…). An empirical tolerance variable is used allowing errors of 1 order per 4 orders (scaling per root)

pasquale90 force-pushed the coef2roots branch from 06df802 to 50f60fe Compare June 13, 2026 10:15

apply micro-opt by replacing divisions with multiplications a) on the…

f4bb673

… loop over subdiabonal elements and b) by avoiding double call of sqrt + replacing division with 2.0 by multipying with 0.5 instead on the extractRoots method

skip updating of the filterstate in case that computed poles are viol…

445fdd6

…ating causality constraint (at least one pole outside the unit circle - magnitude >= 1.0)

pasquale90 closed this Jun 14, 2026

pasquale90 reopened this Jun 14, 2026

	size_t shift_idx {degree}, iter {0};
	size_t shift_idx {degree}, iter {0};
	std::vector<double> v(shift_idx);

Conversation

pasquale90 commented Jun 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

justonem0reuser left a comment

Choose a reason for hiding this comment

Uh oh!

justonem0reuser Jun 10, 2026

Choose a reason for hiding this comment

Uh oh!

pasquale90 Jun 14, 2026

Choose a reason for hiding this comment

Uh oh!

justonem0reuser Jun 10, 2026

Choose a reason for hiding this comment

Uh oh!

justonem0reuser Jun 10, 2026

Choose a reason for hiding this comment

Uh oh!

pasquale90 Jun 14, 2026

Choose a reason for hiding this comment

Uh oh!

justonem0reuser Jun 10, 2026

Choose a reason for hiding this comment

Uh oh!

pasquale90 Jun 14, 2026

Choose a reason for hiding this comment

Uh oh!

ry-gatoni commented Jun 12, 2026

Uh oh!

pasquale90 commented Jun 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

regarding tests

Uh oh!

pasquale90 commented Jun 14, 2026

Uh oh!

pasquale90 commented Jun 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ry-gatoni commented Jun 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

pasquale90 commented Jun 10, 2026 •

edited

Loading

pasquale90 commented Jun 14, 2026 •

edited

Loading

pasquale90 commented Jun 14, 2026 •

edited

Loading